Multicritical Points and Crossover Mediating the Strong Violation of Universality: 
Wang-Landau Determinations in the Random-Bond d = 2 Blume-Capel model 
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The effects of bond randomness on the phase diagram and critical behavior of the square lat- 
tice ferromagnetic Blume-Capel model are discussed. The system is studied in both the pure and 
disordered versions by the same efficient two-stage Wang-Landau method for many values of the 
crystal field, restricted here in the second-order phase transition regime of the pure model. For the 
random-bond version several disorder strengths are considered. We present phase diagram points of 
both pure and random versions and for a particular disorder strength we locate the emergence of the 
enhancement of ferromagnetic order observed in an earlier study in the ex-first-order regime. The 
critical properties of the pure model are contrasted and compared to those of the random model. 
Accepting, for the weak random version, the assumption of the double logarithmic scenario for the 
specific heat we attempt to estimate the range of universality between the pure and random-bond 
models. The behavior of the strong disorder regime is also discussed and a rather complex and yet 
not fully understood behavior is observed. It is pointed out that this complexity is related to the 
ground-state structure of the random-bond version. 

PACS numbers: 75.10.Nr, 05.50.+q, 64.60.Cn, 75.10.Hk 



I. INTRODUCTION 

The effect of quenched randomness on the equilib- 
rium and dynamic properties of macroscopic systems 
is a subject of great theoretical and practical interest. 
It is well known that quenched bond randomness may 
produce drastic changes on phase transitions depending 
on the type of the transition [H-Q. Thus, symmetry- 
breaking first-order transitions are converted to second- 
order phase transitions by infinitesimal bond randomness 
for spatial dimensionality d = 2 0, HJ and by bond ran- 
domness beyond a threshold strength in d > 2 |^|, as 
indicated by general arguments J5[ and in some cases 
by rigorous mathematical work 3j. In particular, this 
rounding effect of first-order transitions has now been 
rigorously established in a unified way in low dimensions 
(d < 2) including a large variety of types of randomness 
in classical and quantum spin systems Q. 

Historically, the effects of disorder on phase transitions 
have been studied in two extreme cases, i.e. in the lim- 
its of weak and strong (near the percolation point) dis- 
order. The first important conjecture, known today as 
the Harris criterion relates the value of the specific 
heat exponent a in a continuous transition with the ex- 
pected effects of uncorrelated weak disorder in ferromag- 
nets. According to the Harris criterion, for continuous 
phase transitions with a negative exponent a, the intro- 
duction of weak randomness is expected to be an irrel- 
evant field and the disordered system to remain in the 
same universality class. On the other hand, the weakly 
disordered system is expected to be in a different univer- 
sality class in the case of a pure system having a positive 
exponent a. Pure systems with a zero specific heat expo- 



nent (a = 0) are marginal cases of Harris criterion (since 
the criterion does not give any information) and their 
study, upon the introduction of disorder, has been of par- 
ticular interest. The paradigmatic model of the marginal 
case is, of course, the general random 2d Ising model 
(random-site, random-bond, and bond-diluted) and this 
model has been extensively investigated and debated 0- 
l24j ]. Several recent studies, both analytical (renormal- 
ization group and conformal field theories) and numeri- 
cal (mainly Monte Carlo (MC) simulations) devoted to 
this model, have provided very strong evidence in favor 
of the so-called logarithmic corrections's scenario. Ac- 
cording to this, the effect of infinitesimal disorder gives 
rise to a marginal irrelevance of randomness and besides 
logarithmic corrections, the critical exponents maintain 
their 2d Ising values. In particular, the specific heat is 
expected to slowly diverge with a double-logarithmic de- 
pendence [Il] - fl4j |. Here, we should mention that there is 
not full agreement in the literature and a different sce- 
nario predicts a negative specific heat exponent a leading 
to a saturating behavior [16| . with a corresponding cor- 
relation length exponent v > 2/d [251 ] . 

In general, a unitary and rigorous physical description 
of critical phenomena in disordered systems still lacks 
and certainly, lacking such a description, the study of 
further models for which there is a general agreement 
in the the behavior of the corresponding pure cases is 
very important. Historically, such a suitable candidate 
for testing the above predictions, that has been also 
quite extensively studied, is the general 2d g-state Potts 
model [l7|, |26I430| | . This model includes the Ising model 
(q = 2), cases of a pure system having continuous transi- 
tions with a positive exponent a (q = 3, 4) and also the 
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large q-cases (q > 4) for which one could observe and try 
to classify the above mentioned softening of first-order 
transitions in 2d models. Another similarly interesting 
candidate, not yet as much studied in the random-bond 
version, is the 2d Blume-Capel (BC) model [3lll32l|. We 
may note here that most of the existing literature on the 
BC model with randomness concerns randomness applied 
to the crystal field and/or spin glass exchange interac- 
tions [33l436| . As it is well known, the pure version of the 
BC model undergoes an Ising-like continuous phase tran- 
sition to an ordered ferromagnetic phase as the tempera- 
ture is lowered for crystal-field couplings less than a tri- 
critical value and a first-order transition for larger values 
of the crystal-field coupling. Therefore, this model pro- 
vides also the opportunity to study two different and very 
interesting topics of the above described effects of disor- 
der in critical phenomena, namely the double-logarithmic 
scenario for the specific heat in the regime where the 2d 
BC model is in the same universality class with the 2d 
random-bond Ising model and also the softening of the 
transition in the first-order regime. Recently the present 
authors [37| have considered this model and provided 
strong numerical evidence clarifying two of the above 
mentioned effects induced in 2d systems by bond random- 
ness. By implementing a two-stage Wang-Landau (WL) 
approach [37H43T ] , we presented essentially exact informa- 
tion on the 2d BC model under quenched bond random- 
ness. In this investigation, we found dramatically dif- 
ferent critical behaviors of the second-order phase transi- 
tions emerging from the first- and second-order regimes of 
the pure BC model and since, these second-order transi- 
tions were found to have different critical exponents, our 
study indicated an interesting strong violation of univer- 
sality [37]. Namely, different sets of critical exponents 
on two segments of the same critical line appeared to 
describe the two regimes: still-second-order and ex-first- 
order. 

In this paper, we extend our earlier work [37| . by im- 
plementing essentially the same two-stage WL approach 
(Sec. HU) and try to give a more complete picture by con- 
centrating in the weak (still-second-order) regime and 
simulate the model for several disorder strengths and 
many values of the crystal-field coupling. The above 
statement means that, effectively we will restrict our 
study to moderate values of the crystal field and moder- 
ate values of the disorder, intending to observe the fron- 
tier between the weak- and the strong-disorder universal- 
ity classes from the disappearance of the expected 2d ran- 
dom Ising universality class behavior. Thus, in Sec. IIIII 
we will produce phase diagram points for the random- 
bond model but also for the pure model, reporting for the 
pure case a comparison with existing estimates in the lit- 
erature. More generally, in carrying out this project we 
have also considered the pure 2d BC model for several 
values of the crystal- field, in the second-order regime, 
observing its finite-size scaling (FSS) behavior. Sec. IIVI 
presents such a comparative study between random and 
pure models concerning the behavior of all thermody- 



namic parameters used in the traditional FSS analysis of 
MC data. This study enables us to observe some pecu- 
liarities of the pure model, due to the onset of tricritical- 
ity, and compare them with the corresponding behavior 
of the random model. Furthermore, we try to focus, un- 
derstand and shed light to the extent of universality of 
the random-bond 2d BC model with the corresponding 
random-bond Ising model, for which the scenario of log- 
arithmic corrections seems to be the strongest option in 
the current literature [24[ • The prediction of the range of 
such universality is far from trivial and the two regimes 
(weak and strong) have many dissimilarities which are 
also reflected in the ground-state structure, as further 
discussed below. The attempt to estimate the range of 
the above mentioned universality is accomplished by the 
novel idea which assumes the truth of the double logarith- 
mic scenario for the specific heat in a suitable restricted 
range. This is presented in Sec. |V] together with some fur- 
ther crucial observations concerning the behavior of the 
strong disorder regime, i.e. the regime where the Ising 
class universality does not apply and the system has a 
rather complex and yet not understood behavior. Our 
conclusions are summarized in Sec. I VI I 



II. DEFINITION OF THE MODELS AND THE 
TWO-STAGE WANG-LANDAU APPROACH 



The pure and random-bond Blume-Capel 
models 



The (pure) BC model [3l|, 133 is defined by the Hamil- 
tonian 



-J 



<ij> 



SiSj 



(1) 



where the spin variables Si take on the values —1,0, 
or +1, < ij > indicates summation over all nearest- 
neighbor pairs of sites, and J > is the ferromag- 
netic exchange interaction. The parameter A is known 
as the crystal-field coupling and to fix the tempera- 
ture scale we set J = 1 and ks = 1- As it is well 
known, this model has been analyzed, besides the orig- 
inal mean- field theory [U, [32[, by a variety of approx- 
imations and numerical approaches. These include the 
real space renormalization group, MC simulations, and 
MC renormalization-group calculations 44|, e-expansion 
renormalization groups [451 ] . high- and low-temperature 
series calculations [46|], a phenomenological FSS analy- 
sis using a strip geometry Ell , and, finally, a recent 
two-parameter WL sampling in rather small lattices of 
linear sizes L < 16 49]. As mentioned already in the in- 
troduction the phase diagram of the model consists of a 
segment of continuous Ising-like transitions at high tem- 
peratures and low values of the crystal field which ends 
at a tricritical point, where it is joined with a second 
segment of first-order transitions between (A t ,T t ) and 
(A = 2,T = 0). 
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The model given by Eq. ([TJ is studied here on the 
square lattice and will be referred to as the pure BC 
model. However, our main focus, on the other hand, is 
the case with bond disorder given by the bimodal distri- 
bution 



P(Jij) 



1 

2 
Ji 



+ J 2 



+ 5(Jij - J 2 )] ; 
Ji > J 2 > ; 



(2) 



■h 
Ji 



so that r reflects the strength of the bond randomness. 
The resulting quenched disordered (random-bond) ver- 
sion of the Hamiltonian defined in Eq. (TTJ) reads now as 



(3) 



A first comparative study between the two versions 
(pure and random) of the 2d BC model has been already 
presented in our earlier paper [37| for two values of the 
crystal-field coupling corresponding to the second-order 
(A = 1) and first-order (A = 1.975) regimes of the pure 
model where for the random version the disorder strength 
r = 0.75/1.25 = 0.6 was chosen in both cases. In the 
next Sections our study on the 2d BC model is extended 
to several values of the crystal field and disorder (listed 
in the table of Sec. IIII[) . Details of our simulations are 
summarized in the next Section together with an up to 
date brief sketch of our entropic scheme. 



B. An outline of our implementation of the 
Wang-Landau approach 

In the last few years we have used an entropic sam- 
pling implementation of the WL algorithm [H, [39| to 
study some simple (iM l4lj , but also some more complex 
systems [1, H3, EH . One basic ingredient of this imple- 
mentation is a suitable restriction of the energy subspace 
for the implementation of the WL algorithm. This was 
originally termed as the critical minimum energy sub- 
space (CrMES) restriction [40j, |41| and it can be carried 
out in many alternative ways, the simplest being that of 
observing the finite-size behavior of the tails of the en ergy 
probability density function (e-pdf) of the system [4lj . 
Complications that may arise in random systems can 
be easily accounted for by various simple modifications 
that take into account possible oscillations in the e-pdf 
and expected sample-to-sample fluctuations of individ- 
ual disorder realizations. In our recent papers [8l. l37l l43j . 
we have presented details of various sophisticated routes 
for the identification of the appropriate energy subspace 
{Ei^Ez) for the entropic sampling of each random real- 
ization. In estimating the appropriate subspace from a 
chosen pseudocritical temperature one should be careful 
to account for the shift behavior of other important pseu- 
docritical temperatures and extend the subspace appro- 
priately from both low- and high-energy sides in order to 



achieve an accurate estimation of all finite-size anoma- 
lies. Of course, taking the union of the corresponding 
subspaces, insures accuracy for the temperature region 
of all studied pseudocritical temperatures. 

The up to date version of our implementation uses a 
combination of several stages of the WL process. First, 
we carry out a starting (or preliminary) multi-range 
(multi-R) stage, in a very wide energy subspace. This 
preliminary stage is performed up to a certain level of 
the WL random walk. The WL refinement is G(E) — >• 
/ * G(E), where G(E) is the density of states (DOS) 
and we follow the usual modification factor adjustment 
fj+i = Wfj and fi — e. The preliminary stage may con- 
sist of the levels : j = 1, . . . ,j = 18 and to improve accu- 
racy the process may be repeated several times. However, 
in repeating the preliminary process and in order to be 
efficient, we use only the levels j — 13, . . . , 18 after the 
first attempt, using as starting DOS the one obtained in 
the first random walk at the level j = 12. From our expe- 
rience, this practice is almost equivalent of simulating the 
same number of independent WL random walks. Also in 
our recent studies we have found out that is much more 
efficient and accurate to loosen up the originally applied 
very strict flatness criteria 53E|- Thus, a variable flat- 
ness process starting at the first levels with a very loose 
flatness criteria and assuming at the level j = 18 the orig- 
inal strict flatness criteria is now days used. After the 
above described preliminary multi-R stage, in the wide 
energy subspace, one can proceed in a safe identification 
of the appropriate energy su bsp ace using one or more al- 
ternatives outlined in Refs. |4(J,|41|. In random systems, 
where one needs to simulate many disorder realizations, 
it is also possible and advisable to avoid the identifica- 
tion of the appropriate energy subspace separately for 
each disorder realization by extrapolating from smaller 
lattices and/or by prediction from preliminary runs on 
small numbers of disorder realizations. In any case, the 
appropriate subspaces should be defined with sufficient 
tolerances. In our implementation we use such advance 
information to proceed in the next stages of the entropic 
sampling. 

The process continues in two further stages (two-stage 
process), using now mainly high iteration levels, where 
the modification factor is very close to unity and there 
is not any significant violation of the detailed balance 
condition during the WL process. These two stages are 
suitable for the accumulation of histogram data (for in- 
stance energy-magnetization histograms), which can be 
used for an accurate entropic calculation of non-thermal 
thermodynamic parameters, such us the order parame- 
ter and its susceptibility [41| . In the first (high-level) 
stage, we follow again a repeated several times (typically 
~ 5 — 10) multi-R WL approach, carried out now only 
in the restricted energy subspace. The WL levels may 
be now chosen as j — 18, 19, 20 and as an appropriate 
starting DOS for the corresponding starting level the av- 
erage DOS of the preliminary stage at the starting level 
may be used. Finally, the second (high-level) stage is ap- 
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plied in the refinement WL levels j = ji, . . . , ji + 3 (typ- 
ically ji — 21), where we usually test both an one-range 
(one-R) or a multi-R approach with large energy inter- 
vals. In the case of the one-R approach we have found 
very convenient and in most cases more accurate to fol- 
low the Belardinelli and Pereyra [Ho| adjustment of the 
WL modification factor according to the rule In / ~ t -1 . 
Finally, it should be also noted that by applying in our 
scheme a separate accumulation of histogram data in the 
starting multi-R stage (in the wide energy subspace) of- 
fers the opportunity to inspect the behavior of all basic 
thermodynamic functions in an also wide temperature 
range and not only in the neighborhood of the finite-size 
anomalies. The approximation outside the dominant en- 
ergy subspace is not of the same accuracy with that of the 
restricted dominant energy subspace but is good enough 
for the observation of the general behavior and provides 
also a route of inspecting the degree of approximation. 

The above described numerical approach was used to 
estimate the properties of a large number of 100 bond 
disorder realizations, for lattice sizes L = 20 — 100 for all 
crystal fields and disorder strengths used in this paper, 
with the exception of the case A — 1.5 and r = 0.5/1.5, 
where 500 disorder realizations were simulated for lat- 
tice sizes L = 20 — 140. For reference and contrast, the 
pure-system's properties were also obtained by the same 
implementation, simulating in each case at least 30 in- 
dependent runs. We close this outline of our numerical 
scheme with some comments concerning statistical er- 
rors and disorder averaging. Even for the larger lattice 
size studied here (L = 100), and depending on the ther- 
modynamic parameter, the statistical errors of the WL 
method were found to be of reasonable magnitude and 
in some cases to be of the order of the symbol sizes, or 
even smaller. This was true for both the pure version and 
the individual random-bond realizations. However, since 
the nature of the present study is qualitative, not aiming 
to an accurate exponent estimation, the WL errors will 
not be presented in our figure illustrations. We also note 
that, for the random-bond version mainly the averages 
over the disorder realizations, denoted as [. . .] av , will be 
considered in the text and their finite-size anomalies, de- 
noted as [. ..]„„, will be used in our FSS attempts. Due 
to very large sample-to-sample fluctuations, mean values 
of individual maxima ([...*]„„) have not been used in 
this study except for illustrative purposes, as in the case 
A = 1.5, r = 0.5/1.5 presented in our last Section. 



III. PHASE DIAGRAMS: PURE AND 
RANDOM-BOND 2D BLUME-CAPEL MODELS 

This Section presents, as mentioned in the introduc- 
tion, phase diagram points for the pure and random-bond 
models and in the case of the pure model compare the 
corresponding phase diagram points with the existing es- 
timates in the literature. This gives also the opportunity 
to observe the reliability of our numerical approach. Fol- 
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FIG. 1: (a) Simultaneous fitting of the six pseudocritical tem- 
perature defined in the text for the random-bond version for 
A = 1.5 and r = 0.9/1.1. The fitting range is L = 20 - 100. 
(b) The same as above simultaneous fitting of the six pseud- 
ocritical temperatures now for the pure model at A = 1.95. 
The solid line shows our estimate of the critical temperature 
T c = 0.659 (with the dotted lines indicating its errors bar- 
riers) and the dashed line the estimate T c — 0.650 given in 
Ref. 0. 



lowing the practice of our earlier paper 37]. we estimate 
phase diagram points by fitting our data to the expected 
power- law shift behavior T = Tc + bL" 1 ^ of several pseu- 
docritical temperatures. The traditionally used specific 
heat and magnetic susceptibility peaks, as well as, the 
peaks corresponding to the following logarithmic deriva- 
tives of the powers n — 1, 2, 4 of the order parameter with 
respect to the inverse temperature K = 1/T (5lT |. 



51n(M n ) _ (M n H) 
dK ~~ (M n ) 



-(H), 



(4) 



and the peak corresponding to the absolute order- 
parameter derivative 



d(\M\) 
dK 



= (\M\H)-(\M\){H), 



(5) 



will be implemented for a simultaneous fitting attempt. 

Such simultaneous fitting attempts are presented in 
Fig. [1] In particular Fig. (Ua) presents the shift be- 
havior for the random-bond 2d BC model in the weak 
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FIG. 2: Phase diagrams of the pure 2d BC model and its 
random version at the disorder strength r = 0.75/1.25. The 
inset illustrates the crossing of the phase boundaries including 
an approximate estimate of the crossing point. 



disorder case r — 0.9/1.1 at the value A = 1.5 and 
Fig. [TJb) illustrates the simultaneous fitting attempted 
for the pure model at A = 1.95, further discussed at 
the end of Sec. IIVBI As noted already, the data fitted 
for the random version of the model are only those of 
the pseudocritical temperatures of the peaks of the aver- 
aged over disorder thermodynamic parameters, as indi- 
cated in the figure, where the asterisk denotes the peak of 
the averaged over disorder parameter [. . .]*„. The alter- 
native route of using averages of individual sample pa- 
rameters gives almost identical estimates. However, in 
cases of strong lack of self-averaging, very large sample- 
to-sample fluctuations may be present. For simplicity, 
in all fitting attempts, the whole range L — 20 — 100 
has been used. Following this practice for the pure 
and random version of the 2d BC several phase dia- 
gram points were produced. Figure [5] shows the result- 
ing phase diagrams using our phase diagram points for 
the pure 2d BC model and its random-bond version for 
the disorder strength r = 0.75/1.25 = 0.6. The tricrit- 
ical point shown is taken from the estimate given by 
Beale (A t , T t )=(1.9655(10), 0.610(5)) [48]. For the dis- 
order strength r = 0.75/1.25 the points of the phase dia- 
gram were chosen with the intension to be able to approx- 
imately locate the emergence of the enhancement of fer- 
romagnetic order observed in our earlier study (37j in the 
ex-first-order regime at A = 1.975, where we found a con- 
siderable increase of the critical temperature by ~ 9%. 
A microscopic explanation of this phenomenon, based 
on the preference Si = ±1 states to the strong-coupling 
connectivity sites has been given in Ref. [37] (see also 
Ref. [Hj]), as a microsegregation process, due to quenched 
bond randomness, that evolves continuously within the 
ferromagnetic and paramagnetic phases. The inset of 
Fig. [3] shows that phase diagram of the random version 
crosses the diagram of the pure model at approximately 



the crossing point (A cross , T cross ) = (1.66(2), 1.02(3)), 
well before the tricritical point. The approximation of 
this point was obtained by interpolation using the four 
points of the phase diagram in the range A = 1.0 — 1.9. 

Table |T] summarizes the phase diagram points ob- 
tained in this paper, by the above described traditional 
FFS method for both the pure and random version of 
the 2d BC model for disorder strengths r = 0.9/1.1, 
r = 0.75/1.25, and r = 0.6/1.4, together with corre- 
sponding phase diagram points given by Beale (48j for 
the pure model. Note that, in Silva et al. [4!| one can 
find an analogous table for the pure 2d BC model includ- 
ing the phase diagram points given by Beale [H| and 
the points produced in their two-parametric WL sam- 
pling [49j]. It can be seen that there is an excellent co- 
incidence of our points with those of Beale [48]. The 
points of Beale |48( are based on the very accurate phe - 
nomenological FSS scheme using a strip geometry [47j, 
whereas our points are obtained via the present simul- 
taneous FSS analysis, based on lattices with linear sizes 
L = 20 — 100. The points of Silva et al. [4!| are based in 
much smaller lattices of the two-parameter WL sampling 
(linear sizes L < 16). However, in the case A = 1.95 our 
phase diagram point does not agree, within errors, with 
that of Beale and for this reason our estimation with a 
rather generous error bars (shown also on the panel) has 
been illustrated in Fig. Qlb). The rest of Table Q] con- 
tains our estimates for the random-bond versions. One 
may note from this table that for small values of A for 
instance A = 1.5 the corresponding critical temperatures 
decrease as the disorder becomes stronger (compare the 
three cases of disorder: r = 0.9/1.1, r = 0.75/1.25, and 
r = 0.6/1.4 at this value of the crystal field). On the 
other hand, for A = 1.9 the trend is reversed as can 
be seen by comparing the corresponding three cases of 
disorder. Apparently this is a kind of reflection of the 
phenomenon of the enhancement of ferromagnetic order 
which appears to influence the geometry of the critical 
surface for the 2d random-bond BC model. 



IV. PHASE TRANSITIONS OF THE PURE AND 
RANDOM-BOND 2D BLUME-CAPEL MODELS 

A. Strong violation of universality: The 
ex-first-order regime 

As pointed out in the introduction, in our recent inves- 
tigation of the 2d random-bond BC model [37[ we found 
dramatically different critical behaviors of the second- 
order phase transitions emerging from the first- and 
second-order regimes of the pure model. Namely, dif- 
ferent sets of critical exponents on two segments of the 
same critical line appeared to describe the two regimes: 
the still-second-order and ex-first-order regimes. The 
study in Ref. [37| was carried out for two values of 
the crystal-field coupling corresponding to the second- 
order (A = 1) and first-order (A = 1.975) regimes of 
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TABLE I: Transition temperatures of the pure and random-bond 2d BC model obtained in this paper. Second column from 
reference [jjfl. third and last entries of third and fifth columns from Ref. [37| . 
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the pure model and for the random version the disor- 
der strength r = 0.75/1.25 was chosen in both cases. 
The strong violation of universality observed appeared 
to be the result of the softening of the first-order tran- 
sition due to bond-randomness. Specifically, it was con- 
cluded that the new strong-disorder universality class is 
well described by a correlation length exponent in the 
range v = 1.30(6) — 1.35(5), and exponent ratios j/f 
and j3/v very close to the Ising values 1.75 and 0.125, re- 
spectively [13| . The above weak universality [HI, HH 
seems to be valid between the Ising-like continuous tran- 
sitions of the random-bond 2d BC model for small values 
of A (A = 1.0 in Ref. 37]) and continuous transitions 
belonging to the strong-disorder universality class. 

Therefore, the strong-disorder universality class may 
be characterized by the above distinct value of the corre- 
lation length exponent and a strong saturation of the spe- 
cific heat. Qualitatively this saturating behavior is quite 
instructive and for illustrative reasons is reproduced here 
in Fig- [3] This figure contrasts, at A — 1.975, the specific 
heat's finite-size behavior of the pure 2d BC model (first- 
order regime) and two disordered cases corresponding to 
disorder strengths r = 0.85/1.15 and r = 0.75/1.25. The 
saturation of the specific heat is very clear in both cases 
of the disorder strength. 

It is of interest to point out here that these findings for 
the strong-disorder universality class appear to be fully 
compatible with the classification of phase transitions in 
disordered systems proposed recently by Wu (HHJ. Ac- 
cording to this classification the strong-disorder transi- 
tion is expected to be inhomogeneous and percolative 
with an expected exponent of the order v = 1.34 
Furthermore, it has been suggested to us by Wu 
that the strong lack of self-averaging of this transition 
stems from the above properties. This violation of self- 
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averaging, together with the strong finite-size effects, 
make the systematic MC approach of the strong-disorder 
regime very demanding, if not impractical. On the other 
hand, the weak regime (or Ising universality regime) suf- 
fers a much weaker lack of self-averaging, by at least a 
factor of ~ 12 [37], and a smooth behavior is observed 
at moderate lattice sizes. Thus, aiming here to observe, 
even approximately, the extent of the involved universal- 
ity classes we carried out our study at moderate values 
of the crystal field and disorder, and found a behavior 
quite convincing from which the frontier of the strong 
universality class can be estimated by observing the dis- 
appearance of the expected 2d random Ising universality 
class. 



B. Pure and random-bond 2d BC model: Range of 
universality with the 2d Ising model 

Let us now proceed with the analysis of our numerical 
data for the disorder strengths and crystal fields given 
in Table U and observe and contrast their FSS behav- 
ior with that of the pure model. Starting this com- 
parative study with the FSS of the specific heat max- 
ima (using for the random-bond version at the strength 
r = 0.75/1.25 the corresponding quantity averaged over 
disorder, i.e. [C]*^), we present in Fig. 0] fitting at- 
tempts for the same range of A for the pure model 
[Fig. Ufa)] and the random-bond version [Fig. HJb)]. As 
indicated by the scales in the x-axis and the functions 
in the corresponding panels, the expected Ising logarith- 
mic divergence has been assumed for the pure model 
C* = C\ + C^lni, whereas the double- logarithmic di- 
vergence [C]* v = C\ + Ci ln(lnL) has been assumed for 
the random version. Although, it is very difficult to ir- 
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FIG. 3: Behavior of the random-bond 2d BC model at 
A = 1.975 from Ref. [37]] • Illustration of the divergence of 
the specific heat of the pure model (first-order regime) and 
the clear saturation of the specific heat for the random-bond 
(open symbols) 2d BC model for two disorder strengths in a 
log-log scale. 



refutably distinguish between a double-logarithmic diver- 
gence and a very weak power-law divergence, the theoret- 
ically well-grounded double-logarithmic scenario applies 
very well [37| and this fact can be observed also now for 
more values of A in Fig. IDJb). There are some further 
features one can observe from Fig. 01 First, for the pure 
model, and in the range A > 1.9, we observe a sud- 
den change in the behavior of the specific heat peaks as 
we approach the tricritical point which is apparently a 
strong cross-over effect. For the random version and the 
same values of A no such strong effects are noticeable 
and most probably the general softening effects of bond 
randomness extends also to the expected cross-over phe- 
nomena between the two different universality classes of 
the random 2d BC model. From Fig. Q|b) the slopes of 
the double-logarithmic fittings appear to obey a rather 
sensible decreasing tendency from which we try in the 
next Section to locate the frontier of strong-disorder uni- 
versality class. 

A second interesting comparison follows now in Fig. [SJ 
where again Fig. [S^a) presents the FSS of the suscepti- 
bility maxima for the pure model, whereas Fig. [SJb) cor- 
responds to the random-bond version at r — 0.75/1.25. 
The influence of the exceptionally large fluctuations in 
the order parameter of the pure model, as we approach 
the tricritical point, is now reflected in the effective ex- 
ponent ratio 7/V, estimated by the simple power law 
X* ~ L 1 ^ . These large fluctuations are a known pe- 
culiarity of the pure model near tricriticality 48]. The 
effective value of the exponent is now closer to the ex- 
pected value at the tricritical point 7/1/ = 1.5 Q than 
to the value of the Ising universality class j/u = 1.75. 
For this reason a separate power-law fitting has been ap- 
plied for the value A = 1.95 in Fig. [SJa). Note here that 
in Fig.[5ja) the simultaneous fitting is applied to the first 
six values of A = — 1.9, whereas in Fig. [5jb) all values 
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FIG. 4: FSS of the specific heat maxima for the pure and 
random-bond 2d BC model at the same values of A. (a) Pure 
model: illustration of linear fittings assuming the expected 
Ising logarithmic divergence, (b) Random model: illustration 
of linear fittings assuming a double logarithmic divergence. 
Note the steady fall of the double logarithmic amplitude Ci. 



of A = — 1.95 are used. However, even so, the com- 
parison of the two panels of Fig. [5j points out the much 
better limiting behavior of the random version towards 
the expected Ising value 7/V = 1.75. Therefore, we may 
convincingly suggest that, for moderate values of crystal 
field and disorder, the Ising universality scenario (with 
possible logarithmic corrections) is well obtained in the 
disordered case. 

Figures [6] and [7] illustrate further alternative routes, 
used commonly in traditional FSS analysis, that pro- 
vide clear evidence to the above suggestion namely that 
the weak-disorder version belong to the 2d Ising class 
for suitable moderate values of disorder and crystal-field 
coupling. Noteworthy is the fact that the estimation of 
the critical exponents via the traditional FSS, such as 
that shown in Figs. [5]- [7] for the random version, yields 
estimates very close to the expected values, i.e. the ex- 
ponents of the 2d Ising model. On the other hand, as 
pointed also earlier, for the pure model, as we approach 
the tricritical point (A = 1.95) the effective exponents 
remind more those expected at the tricritical point than 
those of the 2d Ising model. This is particularly true for 
the correlation length's exponent estimated in Fig. (Ub) 
by the simultaneous fitting of the six pseudocritical tern- 
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FIG. 5: FSS behavior of the susceptibility maxima for the 
pure and random-bond 2d BC model at the same values of 
A. (a) Pure model: simultaneous fitting to a simple power 
law for the first five values of A and a separate fitting close 
to the tricritical point for A = 1.95. (b) Random model: 
simultaneous fitting of the averaged susceptibility peaks for all 
values of A. Note the better behavior for the random model 
and the improved estimation of the exponent ratio ■y/v. 

peratures. As can be seem from this figure the estimate 
for the exponent v, in the case A = 1.95, is closer to the 
expected tricritical value 40/77 = 0.519 •• • [IH, than to 
the 2d Ising value v = 1. 



V. MULTICRITICAL POINTS AND THE 
STRONG DISORDER REGIME 

A. A novel estimation of multicritical points 

From Fig.QJa) one can observe the expected Ising loga- 
rithmic divergence of the specific heat maxima. Avoiding 
the value A = 1.95, which suffers from strong cross-over 
effects, we attempted to estimate the tricritical value of 
the crystal field by fitting the decreasing logarithmic am- 
plitudes C*2 to a suitable power law, as shown in Fig. [3] 
This may appear a naive or questionable idea, since the 
behavior of specific heat data is the Achilles' heel of 
FSS analysis. Yet, Fig. [5] shows that besides the large 
fluctuations (errors) in logarithmic amplitudes C2, one 
could approximately estimate the tricritical crystal field 
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FIG. 6: (a) FSS analysis of the three logarithmic derivatives 
(n = 1,2,4) of the order parameter with respect to tem- 
perature for a particular crystal field A = 0.5 and disorder 
strength r = 0.75/1.25. (b) Traditional FSS analysis of the 
order parameter at the estimated critical temperature for the 
same value of crystal field and disorder strength, as in panel 
(a). 



At f» 1.96(1), as shown in the panel of Fig. |8l 

From Fig. HJb) we observe again that, the specific 
heat maxima, corresponding now to the random-bond 
model at the disorder strength r = 0.75/1.25, are better- 
matched to the double-logarithmic divergence than the 
corresponding specific heat maxima of the pure model 
to a simple logarithmic divergence. Therefore, it appears 
realistic to try to obtain the multicritical point (more pre- 
cisely the value of the crystal field A s d where we expect 
the emergence of strong-disorder regime at a given dis- 
order strength), where the 2d random-bond (BC) Ising 
universality class meets the strong-disorder universality 
class, by fitting the decreasing double-logarithmic am- 
plitudes to a similar power law. Figure [5] presents now 
these fittings and also the estimated multicritical values 
of the crystal field A S( j for the three disorder strengths 
considered in this paper. The statistical errors of the cor- 
responding double logarithmic amplitudes C2 are seen to 
be quite smaller, compared to the corresponding ampli- 
tudes of the pure model, and the estimated values of 
the values of the multicritical field, shown in the panel of 
Fig. [5J appear more convincing. From their general trend 
we observe that as we increase the disorder strength the 
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FIG. 7: Simultaneous fitting to a simple power law of the av- 
eraged peaks corresponding to the absolute order-parameter 
derivative for three values of A and disorder strength r = 0.6. 
The fitting provides an estimate of the exponent (1 — 0)/v 
and as indicated in the panel an alternative estimate for the 
exponent ratio P/u, by assuming v = 1. 



frontier of the strong-disorder universality class moves to 
lower values of the crystal field, namely: (r = 0.9/1.1, 
A sd = 1.963(8)), (r = 0.75/1.25, A sd = 1.955(5)), and 
(r = 0.6/1.4, A sd = 1.879(12)). This behavior seems 
sensible and in our opinion reflects the competition be- 
tween the ferromagnetic random interactions with the 
crystal field, giving a kind of destabilization of the usual 
Ising-like ferromagnetic order. A similar ground-state 
reflection of this competition is a ground-state structure 
of unsaturated ferromagnetic ground states discussed in 
our recent paper [58[. Conversely, as disorder strength is 
decreased, A s d approaches A t , as expected. 



B. Strong disorder regime: The case A = 1.5, 
r = 0.5/1.5 and general observations 

As pointed out earlier, the strong lack of self-averaging, 
together with possible strong finite-size effects, make 
the MC approach to the strong-disorder regime a very 
difficult task. The self-averaging properties along the 
two segments (ex-first-order and still-second-order) of the 
critical line were observed and discussed in Ref. |37J. It 
was shown in this paper that the usual finite-size mea- 
sure [59| of relative variance Rx — Vx/[X^ V , where 
V x = [X 2 ]av - [Xf av (and X = X * is the susceptibil- 
ity maxima), exhibits lack of self-averaging in both cases 
of marginal second-order transition and the transition in 
the strong-disorder or ex-first-order regime. In particu- 
lar, it was shown [37j that the case studied (A = 1.975) of 
the ex-first-order segment gives a much larger effect when 
compared with the still-second-order regime at A = 1 
and the same disorder strength r = 0.75/1.25 by a fac- 
tor of ~ 12. Therefore, phase diagram points very close 
to the frontier of the strong-disorder regime may be the 
worst cases to study, because, besides the very strong 



FIG. 8: An approximate estimation of the tricritical value of 
the crystal field by fitting the decreasing logarithmic ampli- 
tudes Ci of the pure model at suitable values of the crystal- 
field coupling to a suitable power law shown in the figure. 

lack of self-averaging one may also expect large finite-size 
and cross-over effects. However, such cases are elucida- 
tory not only for observing the intrinsic difficulties but 
also for giving us the opportunity to reflect on possible 
links with basic properties of the system as for instance 
the ground-state structure. 

Thus, Fig. [TU1 illustrates two important characteristics 
of the case A = 1.5 and r — 0.5/1.5. In particular, 
Fig. [TUTa) shows the huge sample-to-sample fluctuations 
in all the pseudocritical temperatures used in this pa- 
per. The simulation here was extended to larger lattices 
(L = 20 — 140) and the number of realizations studied 
was 500, to be compared with 100 realizations studied 
in previous Sections (weak- regime) . Besides the enor- 
mous fluctuations, the number of 500 realizations looks 
quite insufficient and strong finite-size effects, reflected 
as expected in the shift behavior of the system, produce 
a completely unsettled behavior. Figure fTUTb 1 ) illustrates 
an unusual deviating behavior between the finite-size be- 
havior of the maxima of the averaged specific heat curves 
([C]*„) and the finite-size behavior of the average of in- 
dividual maxima ([C*] a „), which as shown suffer large 
sample-to-sample fluctuations. The behavior here resem- 
bles in many aspects the well known and still challeng- 
ing specific heat behavior of the 3d random-field Ising 
model [42]. However, a very clear tendency of [C]* av 
for a saturating behavior is observed in Fig. HUf b) and 
we may speculate that this saturating behavior is the 
correct asymptotic behavior for both maxima shown in 
Fig. HOTb). although the behavior of [C*] av will settle 
down only in very large lattice sizes, when the influence 
of the strong finite-size effects on the individual maxima 
will diminish. 

The above illustrations open the possibility that the 
case A = 1.5 and r = 0.5/1.5 is very close to the fron- 
tier of the strong-universality class. This is in accor- 
dance with the trend observed in the previous Section 
and can be further supported by reproducing here some 
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FIG. 9: Estimation of multicritical values of the crystal field, 
where the 2d random-bond (BC) Ising universality class meets 
the strong-disorder universality class. The decreasing double 
logarithmic amplitudes C2 of the random version has been 
fitted to a power law shown and the estimates for the three 
disorder strengths are shown in the panel. 



aspects of the ground-state structure of the 2d random- 
bond BC model. From Fig. [TTJ reproduced here from 
Ref. [58], one can see that approximately at this point 
(A = 1.5, r = 0.5/1.5), the system departs from the ferro- 
magnetic ground state and an unsaturated ground state 
is produced, which is further enhanced with vacant sites 
(sj = 0) as we increase the disorder strength. In the pres- 
ence of bond randomness the competition between the 
ferromagnetic interactions with the crystal field results in 
a destabilization of the ferromagnetic ground state. De- 
pending on the realization, weak clusters exist in T = 
and their points are frozen in the Si — state. This is 
an interesting subject, which is presently under further 
consideration in both 2d and 3d by the present authors. 
The novel behavior illustrated in Fig. [TU] appears now 
as a consequence of the onset of the unsaturated ground 
state at (A = 1.5, r = 0.5/1.5) which is thus related with 
the critical behavior of the 2d random-bond BC model. 
The presented calculation of the ground states has be 
carried out in polynomially bounded computing time by 
mapping the system into a network and searching for a 
minimum cut by using a maximum flow algorithm (see 
for instance Ref. [6(|) and can be easily extended to large 
lattices and also to the 3d BC random-bond model. 



VI. CONCLUSIONS 

By carrying out an extensive two-stage Wang-Landau 
entropic sampling of both the pure and the random- 
bond 2d Blume-Capel model we have produced phase 
diagram points for several disorder strengths. Also for 
the pure model we found an excellent coincidence of 
our points with those of Beale 48]. For a particular 
disorder strength (r = 0.75/1.25) we found that, as a 
result of the enhancement of ferromagnetic order, the 
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FIG. 10: Finite-size behavior for the case A = 1.5 and 
r = 0.5/1.5. (a) Behavior of the averaged pseudocritical 
temperatures, corresponding to individual maxima, with an 
illustration of the huge sample-to-sample fluctuations, (b) 
Behavior of the maxima of the averaged specific heat curves 
([C]at)) and the average of individual maxima ([C*] a u) with 
their large sample-to-sample fluctuations. 



phase diagram of the random version crosses that of the 
pure model at approximately the point (A cross , T cross ) = 
(1.66(2), 1.02(3)), well before the tricritical point. 

The critical properties of the pure model were com- 
pared and contrasted to those of the random model and 
for moderate values of the crystal field and disorder we 
found that the Ising universality scenario (with possible 
logarithmic corrections) is well obtained in the case of the 
random version. Furthermore, accepting in this range 
of couplings the assumption of the double logarithmic 
scenario for the specific heat we estimated multicritical 
points, where the 2d random-bond (BC) Ising universal- 
ity class meets the strong-disorder universality class. 

The behavior of the strong disorder regime was also 
critically discussed. The case A = 1.5 and r = 0.5/1.5 
was extensively studied and based on the observed be- 
havior and on the ground-state observations we suggested 
that this case, most likely, lies on the frontier between the 
weak- and strong-disorder universality classes, suffering 
exceptionally strong finite-size effects and a strong lack 
of self-averaging. 

In conclusion, the present paper investigated some 
difficult and important aspects of the random-bond 2d 
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Blume-Capel model. It has been pointed out that the 
behavior of this system is very interesting and in partic- 
ular the strong-disorder regime may include many further 
challenges and open problems that, at the moment, are 
not fully understood. In our opinion, this is a rather 
complex subject deserving further research. 
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